%% JRD model code
%clear all;
rng('shuffle');

% Specify number of participants to simulate
numParticipants = 14;

% Create a matrix for reference directions and specify num of objects and
% IDs
refDirection = [0,1; 1,0; 0,-1; -1,0];
refDir = 1;
%refDirection = [0,1; 0,-1];
%refDir = 0;

headingConditions = [0,45,90,135,180,225,270,315];
pointConditions = headingConditions;

for p=1:numParticipants
% Parameter: noisy object location
if refDir == 1
    a = optFull(p,1);
    b = optFull(p,2);
else
    a = optHalf(p,1);
    b = optHalf(p,2);
end

l = 1;
q = 20;

% Create a matrix containing the locations of all objects relative to an
% origin (obj id is the index)
locRelOrigin = [-1,0; 0,-1; 0,1; 1,0; 1,1; 1,2; 2,-1; 2,0; 2,1];

% Create cell array; each index is an item containing a n by 2 matrix of 
% vectors representing the distances of each item from that location
% relative to the reference axes. The encoded object locations are then
% computed from those.
for i=1:length(locRelOrigin(:,1))
    for j=1:length(locRelOrigin(:,1))
        locRelObj{i}(j,1) = locRelOrigin(j,1) - locRelOrigin(i,1);
        locRelObj{i}(j,2) = locRelOrigin(j,2) - locRelOrigin(i,2);
    end
end

% Now create an array of test probes (first index = standing location;
% second index = imagined heading; third index = pointing direction) using
% combination algorithm where order matters
n = length(locRelOrigin(:,1)); k = 3;
nk = nchoosek(1:n,k);
testProbe=zeros(0,k);
for i=1:size(nk,1)
    pi = perms(nk(i,:));
    testProbe = unique([testProbe; pi],'rows');
end

% Organize test probes by heading and allocentric pointing direction (round
% to get rid of strange decimals)
for i=1:length(testProbe(:,1))
    if thetaFunc([1,0],locRelObj{testProbe(i,1)}(testProbe(i,2),:)) <= ... 
            thetaFunc([-1,0],locRelObj{testProbe(i,1)}(testProbe(i,2),:))
        testProbe(i,4) = round(thetaFunc([0,1],locRelObj{testProbe(i,1)}(testProbe(i,2),:)));
    else
        testProbe(i,4) = round(180 + thetaFunc([0,-1],locRelObj{testProbe(i,1)}(testProbe(i,2),:)));
    end
    if thetaFunc([1,0],locRelObj{testProbe(i,1)}(testProbe(i,3),:)) <= ... 
            thetaFunc([-1,0],locRelObj{testProbe(i,1)}(testProbe(i,3),:))
        testProbe(i,5) = round(thetaFunc([0,1],locRelObj{testProbe(i,1)}(testProbe(i,3),:)));
    else
        testProbe(i,5) = round(180 + thetaFunc([0,-1],locRelObj{testProbe(i,1)}(testProbe(i,3),:)));
    end
end

% Need to duplicate trials for heading and allocentric condition that only
% appear once given the arrangement of objects
for i=1:length(testProbe(:,1))
    if (testProbe(i,4) == 45 && testProbe(i,5) == 180) || (testProbe(i,4) == ...
            180 && testProbe(i,5) == 45) || (testProbe(i,4) == 45 && ...
            testProbe(i,5) == 270) || (testProbe(i,4) == 270 && testProbe(i,5) ...
            == 45) || (testProbe(i,4) == 135 && testProbe(i,5) == 315) || ...
            (testProbe(i,4) == 315 && testProbe(i,5) == 135)
        testProbe(length(testProbe(:,1)) + 1,:) = testProbe(i,:);
    end
end

    
    % Randomize sample of test probes
    testProbe = testProbe(randperm(size(testProbe,1)),:);
    pointResponse = [];
    
    % Get 128 test probes, a trial for each heading and pointing direction
    count = zeros(8,8);
    newProbeCount = 0;
    for i=1:length(testProbe(:,1))
        for j=1:length(headingConditions)
            if round(testProbe(i,4)) == round(headingConditions(j)) && sum(count(j,:)) < 16
                for k=1:length(pointConditions)
                    if round(testProbe(i,5)) == round(pointConditions(k)) && count(j,k) < 2
                        count(j,k) = count(j,k) + 1;
                        newProbeCount = newProbeCount + 1;
                        newTestProbe(newProbeCount,:) = testProbe(i,:);
                    end
                end
            end
        end
    end
    
    % Randomize new test probes
    newTestProbe = newTestProbe(randperm(size(newTestProbe,1)),:);
            
    % Objects with noisy encoding
    for i=1:length(locRelOrigin(:,1))
        for j=1:length(locRelOrigin(:,1))
            for k=1:length(refDirection(:,1))
                angleHolder(k) = round(thetaFunc(refDirection(k,:),locRelObj{i}(j,:)));           
            end
            [M,I] = min(angleHolder);
            locEncoded{i}(j,1) = locRelObj{i}(j,1) + normrnd(0,(a/(1+b*exp(-(l)*thetaFunc(refDirection(I,:),locRelObj{i}(j,:))))));
            locEncoded{i}(j,2) = locRelObj{i}(j,2) + normrnd(0,(a/(1+b*exp(-(l)*thetaFunc(refDirection(I,:),locRelObj{i}(j,:))))));
            if i > 1
                while i - j > 0    
                    locEncoded{i}(i-j,:) = -(locEncoded{i-j}(i,:));
                    j = j + 1;
                end
            end
        end
    end

    % Now an algorithm to cycle through the test probes and perform the JRD
    % task
    for i=1:length(newTestProbe(:,1))

        % Define standing position, imagined heading, and pointing direction
        standpos = newTestProbe(i,1);
        imagheading = newTestProbe(i,2);
        pointobj = newTestProbe(i,3);

        % Create arrays for heading condition and correct pointing direction
        if thetaFunc([1,0],locRelObj{standpos}(imagheading,:)) <= ... 
                thetaFunc([-1,0],locRelObj{standpos}(imagheading,:))        
            headingCond(i) = thetaFunc([0,1],locRelObj{standpos}(imagheading,:));
        else
            headingCond(i) = 180 + thetaFunc([0,-1],locRelObj{standpos}(imagheading,:));
        end
        % Get rid of strange decimals
        headingCond(i) = round(headingCond(i));
        
        % Get correct pointing angle
        correctPointAngle(i) = atanThetaFunc(locRelObj{standpos}(imagheading,:), ...
            locRelObj{standpos}(pointobj,:));
        % Get rid of strange decimals
        correctPointAngle(i) = round(correctPointAngle(i));

        % Create array for pointing condition
        if thetaFunc([1,0],locRelObj{standpos}(pointobj,:)) <= ... 
                thetaFunc([-1,0],locRelObj{standpos}(pointobj,:))        
            pointingCond(i) = thetaFunc([0,1],locRelObj{standpos}(pointobj,:));
        else
            pointingCond(i) = 180 + thetaFunc([0,-1],locRelObj{standpos}(pointobj,:));
        end
        % Get rid of strange decimals
        pointingCond(i) = round(pointingCond(i));

        % Check if imagined heading is consistent with reference direction
        for j=1:length(refDirection(:,1))
            if round(thetaFunc(locEncoded{standpos}(imagheading,:),refDirection(j,:))) <= q
                pointResponse(i) = atanThetaFunc(refDirection(j,:), ... 
                    locEncoded{standpos}(pointobj,:));
                break
            end
        end

        % If imagined direction is not in line, need to do additional
        % computaions (vecter addition)
        if length(pointResponse) ~= i
            angleVector = locEncoded{standpos}(imagheading,:) + ... 
                locEncoded{imagheading}(pointobj,:);
            pointResponse(i) = atanThetaFunc(locEncoded{standpos}(imagheading,:), ... 
                angleVector);
        end  
        
    end

    % Gather data for pointing errors by heading condition    
    count = zeros(8,1);
    for i=1:length(pointResponse)
        if headingCond(i) == 0
            count(1) = count(1) + 1;
            pointErrorByHeading{1}(count(1)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{1}(count(1)) > 180
                pointErrorByHeading{1}(count(1)) = 360 - pointErrorByHeading{1}(count(1));
            end
        elseif headingCond(i) == 45
            count(2) = count(2) + 1;
            pointErrorByHeading{2}(count(2)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{2}(count(2)) > 180
                pointErrorByHeading{2}(count(2)) = 360 - pointErrorByHeading{2}(count(2));
            end
        elseif headingCond(i) == 90
            count(3) = count(3) + 1;
            pointErrorByHeading{3}(count(3)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{3}(count(3)) > 180
                pointErrorByHeading{3}(count(3)) = 360 - pointErrorByHeading{3}(count(3));
            end
        elseif headingCond(i) == 135
            count(4) = count(4) + 1;
            pointErrorByHeading{4}(count(4)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{4}(count(4)) > 180
                pointErrorByHeading{4}(count(4)) = 360 - pointErrorByHeading{4}(count(4));
            end
        elseif headingCond(i) == 180
            count(5) = count(5) + 1;
            pointErrorByHeading{5}(count(5)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{5}(count(5)) > 180
                pointErrorByHeading{5}(count(5)) = 360 - pointErrorByHeading{5}(count(5));
            end
        elseif headingCond(i) == 225
            count(6) = count(6) + 1;
            pointErrorByHeading{6}(count(6)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{6}(count(6)) > 180
                pointErrorByHeading{6}(count(6)) = 360 - pointErrorByHeading{6}(count(6));
            end
        elseif headingCond(i) == 270
            count(7) = count(7) + 1;
            pointErrorByHeading{7}(count(7)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{7}(count(7)) > 180
                pointErrorByHeading{7}(count(7)) = 360 - pointErrorByHeading{7}(count(7));
            end
        else
            count(8) = count(8) + 1;
            pointErrorByHeading{8}(count(8)) = abs(pointResponse(i) - correctPointAngle(i));
            if pointErrorByHeading{8}(count(8)) > 180
                pointErrorByHeading{8}(count(8)) = 360 - pointErrorByHeading{8}(count(8));
            end
        end
    end
    
    % Gather data for pointing errors by heading direction by pointing
    % condition
    if length(refDirection(:,1)) == 4
        count = zeros(8,2);
        for i=1:length(pointResponse)
            if pointingCond(i) == 0 
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(1,1) = count(1,1) + 1;
                    pointErrorAligned{1}(count(1,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{1}(count(1,1)) > 180
                        pointErrorAligned{1}(count(1,1)) = 360 - pointErrorAligned{1}(count(1,1));
                    end
                else
                    count(1,2) = count(1,2) + 1;
                    pointErrorMisaligned{1}(count(1,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{1}(count(1,2)) > 180
                        pointErrorMisaligned{1}(count(1,2)) = 360 - pointErrorMisaligned{1}(count(1,2));
                    end
                end
            elseif pointingCond(i) == 45
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(2,1) = count(2,1) + 1;
                    pointErrorAligned{2}(count(2,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{2}(count(2,1)) > 180
                        pointErrorAligned{2}(count(2,1)) = 360 - pointErrorAligned{2}(count(2,1));
                    end
                else
                    count(2,2) = count(2,2) + 1;
                    pointErrorMisaligned{2}(count(2,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{2}(count(2,2)) > 180
                        pointErrorMisaligned{2}(count(2,2)) = 360 - pointErrorMisaligned{2}(count(2,2));
                    end
                end
            elseif pointingCond(i) == 90
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(3,1) = count(3,1) + 1;
                    pointErrorAligned{3}(count(3,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{3}(count(3,1)) > 180
                        pointErrorAligned{3}(count(3,1)) = 360 - pointErrorAligned{3}(count(3,1));
                    end
                else
                    count(3,2) = count(3,2) + 1;
                    pointErrorMisaligned{3}(count(3,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{3}(count(3,2)) > 180
                        pointErrorMisaligned{3}(count(3,2)) = 360 - pointErrorMisaligned{3}(count(3,2));
                    end
                end
            elseif pointingCond(i) == 135
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(4,1) = count(4,1) + 1;
                    pointErrorAligned{4}(count(4,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{4}(count(4,1)) > 180
                        pointErrorAligned{4}(count(4,1)) = 360 - pointErrorAligned{4}(count(4,1));
                    end
                else
                    count(4,2) = count(4,2) + 1;
                    pointErrorMisaligned{4}(count(4,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{4}(count(4,2)) > 180
                        pointErrorMisaligned{4}(count(4,2)) = 360 - pointErrorMisaligned{4}(count(4,2));
                    end
                end
            elseif pointingCond(i) == 180
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(5,1) = count(5,1) + 1;
                    pointErrorAligned{5}(count(5,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{5}(count(5,1)) > 180
                        pointErrorAligned{5}(count(5,1)) = 360 - pointErrorAligned{5}(count(5,1));
                    end
                else
                    count(5,2) = count(5,2) + 1;
                    pointErrorMisaligned{5}(count(5,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{5}(count(5,2)) > 180
                        pointErrorMisaligned{5}(count(5,2)) = 360 - pointErrorMisaligned{5}(count(5,2));
                    end
                end
            elseif pointingCond(i) == 225
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(6,1) = count(6,1) + 1;
                    pointErrorAligned{6}(count(6,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{6}(count(6,1)) > 180
                        pointErrorAligned{6}(count(6,1)) = 360 - pointErrorAligned{6}(count(6,1));
                    end
                else
                    count(6,2) = count(6,2) + 1;
                    pointErrorMisaligned{6}(count(6,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{6}(count(6,2)) > 180
                        pointErrorMisaligned{6}(count(6,2)) = 360 - pointErrorMisaligned{6}(count(6,2));
                    end
                end
            elseif pointingCond(i) == 270
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(7,1) = count(7,1) + 1;
                    pointErrorAligned{7}(count(7,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{7}(count(7,1)) > 180
                        pointErrorAligned{7}(count(7,1)) = 360 - pointErrorAligned{7}(count(7,1));
                    end
                else
                    count(7,2) = count(7,2) + 1;
                    pointErrorMisaligned{7}(count(7,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{7}(count(7,2)) > 180
                        pointErrorMisaligned{7}(count(7,2)) = 360 - pointErrorMisaligned{7}(count(7,2));
                    end
                end
            else
                if headingCond(i) == 0 || headingCond(i) == 90 ...
                    || headingCond(i) == 180 || headingCond(i) == 270
                    count(8,1) = count(8,1) + 1;
                    pointErrorAligned{8}(count(8,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorAligned{8}(count(8,1)) > 180
                        pointErrorAligned{8}(count(8,1)) = 360 - pointErrorAligned{8}(count(8,1));
                    end
                else
                    count(8,2) = count(8,2) + 1;
                    pointErrorMisaligned{8}(count(8,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorMisaligned{8}(count(8,2)) > 180
                        pointErrorMisaligned{8}(count(8,2)) = 360 - pointErrorMisaligned{8}(count(8,2));
                    end
                end
            end
        end
    end
    
    if length(refDirection(:,1)) == 2
        count = zeros(8,3);
        for i=1:length(pointResponse)
            if pointingCond(i) == 0
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(1,1) = count(1,1) + 1;
                    pointError0and180{1}(count(1,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{1}(count(1,1)) > 180
                        pointError0and180{1}(count(1,1)) = 360 - pointError0and180{1}(count(1,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(1,2) = count(1,2) + 1;
                    pointError90and270{1}(count(1,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{1}(count(1,2)) > 180
                        pointError90and270{1}(count(1,2)) = 360 - pointError90and270{1}(count(1,2));
                    end
                else
                    count(1,3) = count(1,3) + 1;
                    pointErrorOblique{1}(count(1,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{1}(count(1,3)) > 180
                        pointErrorOblique{1}(count(1,3)) = 360 - pointErrorOblique{1}(count(1,3));
                    end
                end
            end
            if pointingCond(i) == 45
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(2,1) = count(2,1) + 1;
                    pointError0and180{2}(count(2,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{2}(count(2,1)) > 180
                        pointError0and180{2}(count(2,1)) = 360 - pointError0and180{2}(count(2,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(2,2) = count(2,2) + 1;
                    pointError90and270{2}(count(2,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{2}(count(2,2)) > 180
                        pointError90and270{2}(count(2,2)) = 360 - pointError90and270{2}(count(2,2));
                    end
                else
                    count(2,3) = count(2,3) + 1;
                    pointErrorOblique{2}(count(2,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{2}(count(2,3)) > 180
                        pointErrorOblique{2}(count(2,3)) = 360 - pointErrorOblique{2}(count(2,3));
                    end
                end
            end
            if pointingCond(i) == 90
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(3,1) = count(3,1) + 1;
                    pointError0and180{3}(count(3,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{3}(count(3,1)) > 180
                        pointError0and180{3}(count(3,1)) = 360 - pointError0and180{3}(count(3,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(3,2) = count(3,2) + 1;
                    pointError90and270{3}(count(3,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{3}(count(3,2)) > 180
                        pointError90and270{3}(count(3,2)) = 360 - pointError90and270{3}(count(3,2));
                    end
                else
                    count(3,3) = count(3,3) + 1;
                    pointErrorOblique{3}(count(3,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{3}(count(3,3)) > 180
                        pointErrorOblique{3}(count(3,3)) = 360 - pointErrorOblique{3}(count(3,3));
                    end
                end
            end
            if pointingCond(i) == 135
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(4,1) = count(4,1) + 1;
                    pointError0and180{4}(count(4,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{4}(count(4,1)) > 180
                        pointError0and180{4}(count(4,1)) = 360 - pointError0and180{4}(count(4,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(4,2) = count(4,2) + 1;
                    pointError90and270{4}(count(4,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{4}(count(4,2)) > 180
                        pointError90and270{4}(count(4,2)) = 360 - pointError90and270{4}(count(4,2));
                    end
                else
                    count(4,3) = count(4,3) + 1;
                    pointErrorOblique{4}(count(4,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{4}(count(4,3)) > 180
                        pointErrorOblique{4}(count(4,3)) = 360 - pointErrorOblique{4}(count(4,3));
                    end
                end
            end
            if pointingCond(i) == 180
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(5,1) = count(5,1) + 1;
                    pointError0and180{5}(count(5,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{5}(count(5,1)) > 180
                        pointError0and180{5}(count(5,1)) = 360 - pointError0and180{5}(count(5,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(5,2) = count(5,2) + 1;
                    pointError90and270{5}(count(5,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{5}(count(5,2)) > 180
                        pointError90and270{5}(count(5,2)) = 360 - pointError90and270{5}(count(5,2));
                    end
                else
                    count(5,3) = count(5,3) + 1;
                    pointErrorOblique{5}(count(5,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{5}(count(5,3)) > 180
                        pointErrorOblique{5}(count(5,3)) = 360 - pointErrorOblique{5}(count(5,3));
                    end
                end
            end
            if pointingCond(i) == 225
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(6,1) = count(6,1) + 1;
                    pointError0and180{6}(count(6,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{6}(count(6,1)) > 180
                        pointError0and180{6}(count(6,1)) = 360 - pointError0and180{6}(count(6,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(6,2) = count(6,2) + 1;
                    pointError90and270{6}(count(6,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{6}(count(6,2)) > 180
                        pointError90and270{6}(count(6,2)) = 360 - pointError90and270{6}(count(6,2));
                    end
                else
                    count(6,3) = count(6,3) + 1;
                    pointErrorOblique{6}(count(6,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{6}(count(6,3)) > 180
                        pointErrorOblique{6}(count(6,3)) = 360 - pointErrorOblique{6}(count(6,3));
                    end
                end
            end
            if pointingCond(i) == 270
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(7,1) = count(7,1) + 1;
                    pointError0and180{7}(count(7,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{7}(count(7,1)) > 180
                        pointError0and180{7}(count(7,1)) = 360 - pointError0and180{7}(count(7,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(7,2) = count(7,2) + 1;
                    pointError90and270{7}(count(7,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{7}(count(7,2)) > 180
                        pointError90and270{7}(count(7,2)) = 360 - pointError90and270{7}(count(7,2));
                    end
                else
                    count(7,3) = count(7,3) + 1;
                    pointErrorOblique{7}(count(7,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{7}(count(7,3)) > 180
                        pointErrorOblique{7}(count(7,3)) = 360 - pointErrorOblique{7}(count(7,3));
                    end
                end
            end
            if pointingCond(i) == 315
                if headingCond(i) == 0 || headingCond(i) == 180
                    count(8,1) = count(8,1) + 1;
                    pointError0and180{8}(count(8,1)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError0and180{8}(count(8,1)) > 180
                        pointError0and180{8}(count(8,1)) = 360 - pointError0and180{8}(count(8,1));
                    end
                elseif headingCond(i) == 90 || headingCond(i) == 270
                    count(8,2) = count(8,2) + 1;
                    pointError90and270{8}(count(8,2)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointError90and270{8}(count(8,2)) > 180
                        pointError90and270{8}(count(8,2)) = 360 - pointError90and270{8}(count(8,2));
                    end
                else
                    count(8,3) = count(8,3) + 1;
                    pointErrorOblique{8}(count(8,3)) = abs(pointResponse(i) - correctPointAngle(i));
                    if pointErrorOblique{8}(count(8,3)) > 180
                        pointErrorOblique{8}(count(8,3)) = 360 - pointErrorOblique{8}(count(8,3));
                    end
                end
            end
        end
    end
    
    % Get participant averages
    for i=1:length(pointErrorByHeading)
        avgPointErrByHeading(i,p) = mean(pointErrorByHeading{i});
        if length(refDirection(:,1)) == 4
            avgPointErrAligned(i,p) = mean(pointErrorAligned{i});
            avgPointErrMisaligned(i,p) = mean(pointErrorMisaligned{i});
        elseif length(refDirection(:,2)) == 2
            avgPointErr0and180(i,p) = mean(pointError0and180{i});
            avgPointErr90and270(i,p) = mean(pointError90and270{i});
            avgPointErrOblique(i,p) = mean(pointErrorOblique{i});
        end       
    end
end


% Now get averages over participants
for i=1:length(avgPointErrByHeading(:,1))
    avgPointErrHeadingCollapsed(i) = mean(avgPointErrByHeading(i,:));
    if length(refDirection(:,1)) == 4
        avgPointErrAlignedCollapsed(i) = mean(avgPointErrAligned(i,:));
        avgPointErrMisalignedCollapsed(i) = mean(avgPointErrMisaligned(i,:));
    elseif length(refDirection(:,1)) == 2
        avgPointErr0and180Collapsed(i) = mean(avgPointErr0and180(i,:));
        avgPointErr90and270Collapsed(i) = mean(avgPointErr90and270(i,:));
        avgPointErrObliqueCollapsed(i) = mean(avgPointErrOblique(i,:));
    end
end

% Plot results, first by heading condition, then by pointing condition
figure;
plot(headingConditions,avgPointErrHeadingCollapsed, '-s', 'MarkerSize', 8);
xlabel('Heading Condition');
ylabel('Pointing Error');
xlim([-45 360]);
ylim([0 50]);
xticks(0:45:315);
yticks(0:10:50);

if length(refDirection(:,1)) == 4
    figure;
    hold on
    plot(pointConditions,avgPointErrAlignedCollapsed, '-s', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrMisalignedCollapsed, '-s', 'MarkerSize', 8);
    hold off
    xlabel('Pointing Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    legend('Aligned','Misaligned');
elseif length(refDirection(:,1)) == 2
    figure;
    hold on
    plot(pointConditions,avgPointErr0and180Collapsed, '-s', 'MarkerSize', 8);
    plot(pointConditions,avgPointErr90and270Collapsed, '-s', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrObliqueCollapsed, '-s', 'MarkerSize', 8);   
    hold off
    xlabel('Pointing Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    legend('0-180','90-270','45-135-225,315');
end

%% Organize observed data
Exp2matrix = readmatrix('Experiment 2 (LayM).xlsx');
if length(refDirection(:,1)) == 4
    Exp2matrixFull = Exp2matrix(1:(length(Exp2matrix(:,1))/2),:);
    % For easy reading, we'll do heading and allocentric effects separately
    % First, heading
    count = zeros(8,1);
    for i = 1:length(Exp2matrixFull(:,1))
        if Exp2matrixFull(i,5) == 0
            count(1) = count(1) + 1;
            pointErrorByHeadingObs{1}(count(1)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 45
            count(2) = count(2) + 1;
            pointErrorByHeadingObs{2}(count(2)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 90
            count(3) = count(3) + 1;
            pointErrorByHeadingObs{3}(count(3)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 135
            count(4) = count(4) + 1;
            pointErrorByHeadingObs{4}(count(4)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 180
            count(5) = count(5) + 1;
            pointErrorByHeadingObs{5}(count(5)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 225
            count(6) = count(6) + 1;
            pointErrorByHeadingObs{6}(count(6)) = Exp2matrixFull(i,13);
        elseif Exp2matrixFull(i,5) == 270
            count(7) = count(7) + 1;
            pointErrorByHeadingObs{7}(count(7)) = Exp2matrixFull(i,13);
        else
            count(8) = count(8) + 1;
            pointErrorByHeadingObs{8}(count(8)) = Exp2matrixFull(i,13);
        end
    end
    
    % Now, allocentric
    count = zeros(8,2);
    for i=1:length(Exp2matrixFull(:,1))
        if Exp2matrixFull(i,14) == 0 
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(1,1) = count(1,1) + 1;
                pointErrorAlignedObs{1}(count(1,1)) = Exp2matrixFull(i,13);
            else
                count(1,2) = count(1,2) + 1;
                pointErrorMisalignedObs{1}(count(1,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 45
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(2,1) = count(2,1) + 1;
                pointErrorAlignedObs{2}(count(2,1)) = Exp2matrixFull(i,13);
            else
                count(2,2) = count(2,2) + 1;
                pointErrorMisalignedObs{2}(count(2,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 90
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(3,1) = count(3,1) + 1;
                pointErrorAlignedObs{3}(count(3,1)) = Exp2matrixFull(i,13);
            else
                count(3,2) = count(3,2) + 1;
                pointErrorMisalignedObs{3}(count(3,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 135
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(4,1) = count(4,1) + 1;
                pointErrorAlignedObs{4}(count(4,1)) = Exp2matrixFull(i,13);
            else
                count(4,2) = count(4,2) + 1;
                pointErrorMisalignedObs{4}(count(4,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 180
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(5,1) = count(5,1) + 1;
                pointErrorAlignedObs{5}(count(5,1)) = Exp2matrixFull(i,13);
            else
                count(5,2) = count(5,2) + 1;
                pointErrorMisalignedObs{5}(count(5,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 225
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(6,1) = count(6,1) + 1;
                pointErrorAlignedObs{6}(count(6,1)) = Exp2matrixFull(i,13);
            else
                count(6,2) = count(6,2) + 1;
                pointErrorMisalignedObs{6}(count(6,2)) = Exp2matrixFull(i,13);
            end
        elseif Exp2matrixFull(i,14) == 270
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(7,1) = count(7,1) + 1;
                pointErrorAlignedObs{7}(count(7,1)) = Exp2matrixFull(i,13);
            else
                count(7,2) = count(7,2) + 1;
                pointErrorMisalignedObs{7}(count(7,2)) = Exp2matrixFull(i,13);
            end
        else
            if Exp2matrixFull(i,5) == 0 || Exp2matrixFull(i,5) == 90 ...
                || Exp2matrixFull(i,5) == 180 || Exp2matrixFull(i,5) == 270
                count(8,1) = count(8,1) + 1;
                pointErrorAlignedObs{8}(count(8,1)) = Exp2matrixFull(i,13);
            else
                count(8,2) = count(8,2) + 1;
                pointErrorMisalignedObs{8}(count(8,2)) = Exp2matrixFull(i,13);
            end
        end
    end
    
    % Now get averages
    for i=1:length(pointErrorByHeadingObs)
        avgPointErrorByHeadingObs(i) = mean(pointErrorByHeadingObs{i});
        avgPointErrorAlignedObs(i) = mean(pointErrorAlignedObs{i});
        avgPointErrorMisalignedObs(i) = mean(pointErrorMisalignedObs{i});
    end
    
    % Need to put allocentric pointing direction conditions in one vector
    allocentricPointEffectMdl = [avgPointErrAlignedCollapsed avgPointErrMisalignedCollapsed];
    allocentricPointEffectObs = [avgPointErrorAlignedObs avgPointErrorMisalignedObs];

    HeadingEffectFit = corrcoef(avgPointErrorByHeadingObs,avgPointErrHeadingCollapsed)
    AllocentricEffectFit = corrcoef(allocentricPointEffectObs,allocentricPointEffectMdl)
    
    % Plot heading effect
    figure;
    hold on
    plot(headingConditions,avgPointErrHeadingCollapsed, '-sb', 'MarkerSize', 8);
    plot(headingConditions,avgPointErrorByHeadingObs, '--sb', 'MarkerSize', 8);
    hold off
    xlabel('Heading Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    
    % Plot pointing effect
    figure;
    hold on
    plot(pointConditions,avgPointErrAlignedCollapsed, '-sb', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrMisalignedCollapsed, '-sr', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrorAlignedObs, '--sb', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrorMisalignedObs, '--sr', 'MarkerSize', 8);
    hold off
    xlabel('Pointing Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    
else
    Exp2matrixHalf = Exp2matrix((length(Exp2matrix(:,1))/2)+1:length(Exp2matrix(:,1)),:);
    % For easy reading, we'll do heading and allocentric effects separately
    % First, heading
    count = zeros(8,1);
    for i = 1:length(Exp2matrixHalf(:,1))
        if Exp2matrixHalf(i,5) == 0
            count(1) = count(1) + 1;
            pointErrorByHeadingObs{1}(count(1)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 45
            count(2) = count(2) + 1;
            pointErrorByHeadingObs{2}(count(2)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 90
            count(3) = count(3) + 1;
            pointErrorByHeadingObs{3}(count(3)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 135
            count(4) = count(4) + 1;
            pointErrorByHeadingObs{4}(count(4)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 180
            count(5) = count(5) + 1;
            pointErrorByHeadingObs{5}(count(5)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 225
            count(6) = count(6) + 1;
            pointErrorByHeadingObs{6}(count(6)) = Exp2matrixHalf(i,13);
        elseif Exp2matrixHalf(i,5) == 270
            count(7) = count(7) + 1;
            pointErrorByHeadingObs{7}(count(7)) = Exp2matrixHalf(i,13);
        else
            count(8) = count(8) + 1;
            pointErrorByHeadingObs{8}(count(8)) = Exp2matrixHalf(i,13);
        end
    end
    
    % Now, allocentric
    count = zeros(8,3);
    for i=1:length(Exp2matrixHalf(:,1))
        if Exp2matrixHalf(i,14) == 0
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(1,1) = count(1,1) + 1;
                pointError0and180Obs{1}(count(1,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(1,2) = count(1,2) + 1;
                pointError90and270Obs{1}(count(1,2)) = Exp2matrixHalf(i,13);
            else
                count(1,3) = count(1,3) + 1;
                pointErrorObliqueObs{1}(count(1,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 45
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(2,1) = count(2,1) + 1;
                pointError0and180Obs{2}(count(2,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(2,2) = count(2,2) + 1;
                pointError90and270Obs{2}(count(2,2)) = Exp2matrixHalf(i,13);
            else
                count(2,3) = count(2,3) + 1;
                pointErrorObliqueObs{2}(count(2,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 90
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(3,1) = count(3,1) + 1;
                pointError0and180Obs{3}(count(3,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(3,2) = count(3,2) + 1;
                pointError90and270Obs{3}(count(3,2)) = Exp2matrixHalf(i,13);
            else
                count(3,3) = count(3,3) + 1;
                pointErrorObliqueObs{3}(count(3,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 135
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(4,1) = count(4,1) + 1;
                pointError0and180Obs{4}(count(4,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(4,2) = count(4,2) + 1;
                pointError90and270Obs{4}(count(4,2)) = Exp2matrixHalf(i,13);
            else
                count(4,3) = count(4,3) + 1;
                pointErrorObliqueObs{4}(count(4,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 180
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(5,1) = count(5,1) + 1;
                pointError0and180Obs{5}(count(5,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(5,2) = count(5,2) + 1;
                pointError90and270Obs{5}(count(5,2)) = Exp2matrixHalf(i,13);
            else
                count(5,3) = count(5,3) + 1;
                pointErrorObliqueObs{5}(count(5,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 225
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(6,1) = count(6,1) + 1;
                pointError0and180Obs{6}(count(6,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(6,2) = count(6,2) + 1;
                pointError90and270Obs{6}(count(6,2)) = Exp2matrixHalf(i,13);
            else
                count(6,3) = count(6,3) + 1;
                pointErrorObliqueObs{6}(count(6,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 270
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(7,1) = count(7,1) + 1;
                pointError0and180Obs{7}(count(7,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(7,2) = count(7,2) + 1;
                pointError90and270Obs{7}(count(7,2)) = Exp2matrixHalf(i,13);
            else
                count(7,3) = count(7,3) + 1;
                pointErrorObliqueObs{7}(count(7,3)) = Exp2matrixHalf(i,13);
            end
        end
        if Exp2matrixHalf(i,14) == 315
            if Exp2matrixHalf(i,5) == 0 || Exp2matrixHalf(i,5) == 180
                count(8,1) = count(8,1) + 1;
                pointError0and180Obs{8}(count(8,1)) = Exp2matrixHalf(i,13);
            elseif Exp2matrixHalf(i,5) == 90 || Exp2matrixHalf(i,5) == 270
                count(8,2) = count(8,2) + 1;
                pointError90and270Obs{8}(count(8,2)) = Exp2matrixHalf(i,13);
            else
                count(8,3) = count(8,3) + 1;
                pointErrorObliqueObs{8}(count(8,3)) = Exp2matrixHalf(i,13);
            end
        end
    end
    
    % Now get averages
    for i=1:length(pointErrorByHeadingObs)
        avgPointErrorByHeadingObs(i) = mean(pointErrorByHeadingObs{i});
        avgPointError0and180Obs(i) = mean(pointError0and180Obs{i});
        avgPointError90and270Obs(i) = mean(pointError90and270Obs{i});
        avgPointErrorObliqueObs(i) = mean(pointErrorObliqueObs{i});
    end
    
    allocentricPointEffectMdl = [avgPointErr0and180Collapsed avgPointErr90and270Collapsed ... 
        avgPointErrObliqueCollapsed];
    allocentricPointEffectObs = [avgPointError0and180Obs avgPointError90and270Obs ... 
        avgPointErrorObliqueObs];
    
    HeadingEffectFit = corrcoef(avgPointErrorByHeadingObs,avgPointErrHeadingCollapsed)
    AllocentricEffectFit = corrcoef(allocentricPointEffectObs,allocentricPointEffectMdl)
    
    % Plot heading effect
    figure;
    hold on
    plot(headingConditions,avgPointErrHeadingCollapsed, '-sb', 'MarkerSize', 8);
    plot(headingConditions,avgPointErrorByHeadingObs, '--sb', 'MarkerSize', 8);
    hold off
    xlabel('Heading Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    
    % Plot pointing effect
    figure;
    hold on
    plot(pointConditions,avgPointErr0and180Collapsed, '-sb', 'MarkerSize', 8);
    plot(pointConditions,avgPointErr90and270Collapsed, '-sr', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrObliqueCollapsed, '-sy', 'MarkerSize', 8);
    plot(pointConditions,avgPointError0and180Obs, '--sb', 'MarkerSize', 8);
    plot(pointConditions,avgPointError90and270Obs, '--sr', 'MarkerSize', 8);
    plot(pointConditions,avgPointErrorObliqueObs, '--sy', 'MarkerSize', 8);
    hold off
    xlabel('Pointing Condition');
    ylabel('Pointing Error');
    xlim([-45 360]);
    ylim([0 50]);
    xticks(0:45:315);
    yticks(0:10:50);
    
end
